{FLOAT SP_int_x, SM_int_x, ONEoverSPMSM_int_x, SP_int_y, SM_int_y, ONEoverSPMSM_int_y;
FLOAT  V_p_int_x, rN_int_x, V_p_int_y, rN_int_y;

#include"wave_2d_eq_eq.c"


 
 SP_int_x= sqrt(1.0 + x*x)*D;

 rN_int_x = nx*((phi_x - phi_x_n_x) + x*y/(1+x*x)*(phi_y - phi_y_n_x));
 V_p_int_x = 0.5*((phi_t - phi_t_n_x) / sqrt(1.0 + x*x)/D + rN_int_x);

 SP_int_y= sqrt(1.0 + y*y)*D; 

 rN_int_y = ny*((phi_y - phi_y_n_y) + x*y/(1+y*y)*(phi_x - phi_x_n_y));
 V_p_int_y = 0.5*((phi_t - phi_t_n_y) / sqrt(1.0 + y*y)/D + rN_int_y);

	     (*derivs_ptr).u[PHI]GRD_IND = phi_t;
	     (*derivs_ptr).u[PHI_T]GRD_IND = (phi_t_t 
			   - Tau_p_int_x*SP_int_x * sqrt(1.0 + x*x)*D*V_p_int_x
			   - Tau_p_int_y*SP_int_y * sqrt(1.0 + y*y)*D*V_p_int_y
						);
	     (*derivs_ptr).u[PHI_X]GRD_IND = (phi_x_t - Tau_p_int_x*nx*SP_int_x*V_p_int_x);
	     (*derivs_ptr).u[PHI_Y]GRD_IND = (phi_y_t - Tau_p_int_y*ny*SP_int_y*V_p_int_y);
}
